library(leaflet)
library(tidyverse)
library(RSocrata)
library(leaflet)
library(sf)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
In this document we get the parks data from the New York City open data repository and do some cleaning. First, we get the parks dataset. The original data dictionary is available at on Github
parks_url <- "https://data.cityofnewyork.us/api/geospatial/y6ja-fw4f?method=export&format=GeoJSON"
parks <- st_read(parks_url, stringsAsFactors = FALSE) %>%
mutate(shape_area = as.numeric(shape_area))
## Reading layer `OGRGeoJSON' from data source `https://data.cityofnewyork.us/api/geospatial/y6ja-fw4f?method=export&format=GeoJSON' using driver `GeoJSON'
## Simple feature collection with 12491 features and 10 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -74.25537 ymin: 40.49613 xmax: -73.70182 ymax: 40.91138
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
The New York City open data portal also contains a list of designated trails within parks.
trails <- read_csv("https://data.cityofnewyork.us/api/views/vjbm-hsyr/rows.csv?accessType=DOWNLOAD")
## Parsed with column specification:
## cols(
## Park_Name = col_character(),
## Width_ft = col_character(),
## Class = col_character(),
## Surface = col_character(),
## Gen_Topog = col_character(),
## Difficulty = col_double(),
## Date_Collected = col_character(),
## Trail_Name = col_character(),
## ParkID = col_character(),
## TrailMarkersInstalled = col_character(),
## SHAPE = col_character()
## )
trails2 <- jsonlite::fromJSON("https://www.nycgovparks.org/bigapps/DPR_Hiking_001.json")
parks_with_trails <- c(trails$ParkID, trails2$Prop_ID)
The map of Flushing Meadows in Queens below shows some of the difficulties of this dataset: - First, the park is constructed of multiple outlines that limits the size of the total park. This makes Flushing Meadows appear smaller than Central Park, which is a single large park. - Second, the dataset includes tennis parks and baseball fields as separate facilities. We want to include these as attributes of their containing park rather than independent facilities. - Third, parking lots for stadia (in this case Citi Field) are included as open park space. We should remove these as they are categorically unlike other green spaces. - Fourth, important green spaces like cemeteries are in a different dataset.
leaflet(st_buffer(parks %>% filter(grepl("Flushing Meadows", park_name)), dist = 0.000001)) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(label = ~as.character(park_name), stroke = 0.5)
## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
## dist is assumed to be in decimal degrees (arc_degrees).
First, we manually remove the park facilities corresponding to stadia and their associated parking lots, and then tag facilities that are not explicitly park boundaries.
notparks <- c(9498000833, 14498000872, 14498000794, # citi field, Yankee stadium, garage
9498000631) # USTA Billie Jean Tennis Center
parks <- parks %>%
filter(!(source_id %in% notparks)) %>%
mutate(
courts = ifelse(grepl("Courts", landuse), TRUE, FALSE),
playgrounds = ifelse(grepl("Playground", landuse), TRUE, FALSE),
trails = ifelse(parknum %in% parks_with_trails, TRUE, FALSE)
)
With those areas removed, we build a 30-foot buffer around each park and dissolve interior boundaries based on the park number. Tennis courts, etc. become boolean variables indicating whether the park has the facilities or not. This will make parks that are effectively the same facility a single polygon feature. We also re-calculate the area of the parks in acres and remove parks smaller than half an acre.
dissolved <- parks %>%
filter(!is.na(parknum)) %>%
group_by(parknum) %>%
filter(shape_area > 0.5 * 43560 | courts) %>%
summarise(
area = sum(shape_area),
courts = any(courts),
playgrounds = any(playgrounds),
trails = any(trails),
park_name = park_name[1]
) %>%
mutate(sqft = st_area(.), acres = units::set_units(sqft, acres)) %>%
st_transform(3628)
leaflet(st_transform(dissolved, 4326) ) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(stroke = 0.5,)